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SIGNAL PROCESSING 

The present invention relates to the reconstruction of 
signals. There are many applications, such as radar, 
sonar, acoustic data, spectroscopy, geophysical, and image 
signal processing, in which it is desirable to reconstruct 
signals from given data. 

In many of these situations the effect on the signals 
of noise and the particular characteristic of the system 
generating the signal are known or can be approximated 
using appropriate mathematical models. In these 

situations, Bayesian reconstruction methods have often been 
applied to reconstruct the signal from given data. These 
methods can work well. For example one Bayesian 
reconstruction approach, known as the maximum entropy 
method, is known to work well. Usually, the maximum 
entropy method (MEM) , is only be applied to the 
reconstruction of signals that are strictly non-negative, 
and which are not complex. Given this, it is generally not 
possible to change the basis in any vector space 
representing the data during any reconstruction. Being 
unable to do this results in computationally very expensive 
reconstruction processes having to be employed, because of 
the requirement for very large calculations having an 
extremely large number of variables to be determined. 

For example, employment of such an MEM to a stack of 
twenty images from a microscope takes in the region of 

TM TM 

fifty minutes using a standard INTEL Pentium 200MHz 
processor. 

The present invention is directed towards improving 
the reconstruction of signals from given data so that such 
reconstruction can be performed within a reduced time 
frame. 

According to the present invention there is provided 
a method of reconstructing a signal from a given set of 
data, with a prediction function representing a predictable 
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effect on the signal, and a noise function representing 
unpredictable noise , the method comprising the steps of: 

altering the coordinate basis of the data and signal 
from an original coordinate basis in order to produce a 
5 prediction function having a reduced set of variables; 

performing a Bayesian reconstruction capable of 
operation of positive , negative , and complex signal values 
to produce a reconstruction signal; and 

converting the reconstruction signal back into the 
10 original coordinate basis to generate a signal. 

The Bayesian reconstruction may be performed using a 
Fourier basis, or may use a wavelet basis. 

The Bayesian reconstruction may employ the maximum 
entropy method, in which case the method may employ an 
15 evaluation parameter, oc, which may be determined from a 
prior reconstruction, set at a fixed value, or determined 
during the reconstruction step. 

The signal to be reconstructed may be an image signal , 
or may be a radar, sonar, or acoustic data signal. 
20 Alternatively, it may be a signal from spectroscopy or a 
geophysical signal. 

By employing the method of the present invention, an 
example stack of 20 microscope images takes approximately 

TM TM 

45 seconds to reconstruct using a standard INTEL Pentium 
25 200 MH Z processor. 

One example of the present invention will now be 
described. 



* 
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Bayesian reconstruction methods have been applied to 
numerous problems in a wide variety of fields.; In their 
standard form, however, they can be very computationally 
intensive, since they generally require the numerical 
maximisation of a complicated function of many variables. 
For example, in image reconstruction problems, it is not 
unusual for the number of variables to be ' 10 6 . 
Furthermore, one of the most popular Bayesian 
reconstruction algorithms is the maximum-entropy method 
(MEM), which can only by applied to the reconstruction of 
signals that are strictly non-negative (see below). The 
method can, however, be extended to signals that can take 
both positive and negative values. We develop the MEM 
approach so that it can be applied to the reconstruction of 
signals that can take positive, negative or complex values. 
As a result, this enables the use of similarity 
transformations in the reconstruction algorithm so that 
calculations can be performed in an alternative "basis" 
this is more appropriate to the problem under 
consideration. Specifically, the basis is chosen so that 
signal is reconstructed by performing a large number of 
numerical maximisations of low dimensionality, rather than 
a single maximisation of high dimensionality: This results 
in a significant increase in speed. Indeed, in the example 
outlined below, the speed of the reconstruction algorithm 
is increased by a factor of about 100. 

The standard Bayesian reconstruction techniques 
mentioned above are typically applied to a given data set 
d(y) in order to reconstruct some underlying signal s(x) . 
Here, y denotes the space over which the data are defined 
and x denotes the space of the signal, which may in general 
be distinct from x and need not have the same number of 
dimensions. For example, given some data stream d(t), which 
consists of the measured values of some quantity as a 
function of time t, we may wish to reconstruct the two- 
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dimensional spatial variation of some other quantity (or 

signal) s(x,y) . 

For most measurements it is convenient (or necessary) 

to digitise the data and the signal (for example if either 

is to be stored/analyzed on a digital computer). We may 

therefore denote the data by the vector d with N d 

components, where N d is the number of data samples. 

Similarly, we denote the signal by the vector s of length 

Ns . where N s is the number of points at which we wish to 

reconstruct the signal. 

in general, we may express the data vector as some 

function e of the signal vector, i.e. 

d=e(s) 

The function e can be non-linear and specifies the effect 
of the measuring apparatus on the signal that we wish to 
reconstruct. It is customary to divide this function into 
the predictable effect of the measuring apparatus on the 
signal and the stochastic noise part due to inherent 
inaccuracies in the measurement process. In this case, we 

» 

may write 

d = *(s) + <. (1) 

Where * denotes the predictable response of the 
apparatus to the signal and < is a vector of length N d 
that contains any stochastic noise contributions to the 

data . 

The Bayesian approach to reconstructing the signal 
is to calculate the estimator * that maximises the 
posterior probability Pr(sd) . This is given by Bayes' 

theorem as Pr(d|s) Pr(s) 

Pr(s|d) = — p^-fj . 

where Pr(d|s) is the likelihood of obtaining the 
data given the signal, Pr(s) is the prior probability, 
and the evidence Pr(d) can be considered merely as a 
normalisation constant. Thus, in order to obtain the 
(Bayesian estimator of the signal vector, we must 
maximise the product the product Pr(d|s) P(s) of the 
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likelihood function and the prior. A thorough discussion 
of Bayesian analysis techniques is given by Sivia (1996). 

The likelihood function describes the statistics of 
noise contribution t to the data. This function may 
5 take any form appropriate to the noise statistics. It is 
convenient to define the log -likelihood function 
L (s) = In[Pr(d|s)] so that the likelihood function may be 
written as Pr(d|s) = exp[L(s)]. As an example, if the 
noise on the data is Gaussian-distributed and described 
10 by the noise covariance matrix N, then the likelihood 
function takes the form 

Pr(d|s) « exph^N-^] 

cc expHtd-^OrN-Md-*^))]. (2) 



15 



20 



where in the second line we have used (1). In this 
case, the log- likelihood function is simply minus one 
half of the standard x 3 misfit statistic, i.e. 

Us) = 

The prior distribution Pr(s) codifies our knowledge 
of the underlying signal before acquiring the data. If 
we have some advance knowledge of the statistical 
properties of the signal then this is contained in the 
prior. For example, if we assume the signal to be 
described by a Gaussian random field with a covariance 
matrix C, then the prior takes the form 

Pr(s) « exp(-^s T C -1 s). 

3 0 indeed, if the prior is assumed to have this form 

and the likelihood is also Gaussian, as in (2), then the 
Bayesian estimator § is obtained by maximising their 
product is identical to the standard Wiener filter 
solution. An introduction to the Wiener filter technique 

35 is given by Press et al (1994). 



25 



SUBSTITUTE SHEET (RULE 26) 



WO 00/1 4649 PCT/GB99/Q271 5 

6 " — 

It is clear, however, that although the noise 
contribution * to the data may often be Gaussian- 
distributed, the assumption of a Gaussian form for the 
prior is not valid for a general signal. If the joint 
5 probability distribution of the elements of the signal 
vector is known then it should be used as the prior. 
This is almost always impossible, however, and we instead 
investigate the assignment of a prior applicable to 
general signals that is based on information-theoretic 
10 considerations alone. Using very general notions of 
subset independence, coordinate invariance and system 
independence, it may be shown that the prior probability 
Pr(s) should take the form 



15 



25 



30 



35 



p r (s) « exp[cr5(s, m)], ( 3 ) 



Where the dimensional constant a . depends on the 
scaling of the problem and may be considered as a 
regularising parameter, and ro is a model to which the 
20 Bayesian reconstruction defaults in the absence of any 
data and is usually set to a small constant value. The 
function S(s, m) is the cross-entropy of the signal and 
model vectors and is given by 



S(s,m) = £L-m„-* n ln(^)]. (4) 

Where N s , is the length of the (digitised) signal 
vector. A derivation of this result is given by Skilling 
(1988). By combining the entropic prior with the 
likelihood function, the Bayesian estimator of the signal 
is found by maximising with respect to s the posterior 
probability, which now takes the form 

Pr(s|d) « exp[L(s)]exp(aS(s, m)] = exp[L(s) + aS(s,m)]. 

» 

• * 

* 

Thus, maximising this probability distribution is 
equivalent to maximising the function 
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F(s) = L(s) + aS{s,m), (5) 

.. • and this forms the basis of the maximum entropy 
method (MEM) . 

5 

The maximum-entropy method has been applied to a 
wide range of signal reconstruction problems. In its 
standard form, however, it can by very computationally 
intensive. The function F(s) is in general a complicated 

10 function of the components s n of the signal vector and so 
a numerical maximisation of the F(s) must be perform over 
this N a -dimensional space. It is not unusual for N s to 
be of the order N s ~ 10 6 , particularly in image 
reconstruction problems. Moreover, the standard MEM 

15 approach is only applicable to signals that are strictly 
non-negative, as is clear from the presence of the 
logarithmic term in the expression (4) for the entropy. 

Nevertheless, it is possible to extend the MEM 
so that it can be applied to the reconstruction of 

20 signals that can take both positive and negative values. 
The definition of the entropy for positive/negative 
signals with certain special properties was first 
presented by Gull & Skilling (1990). The generalisation 
to arbitrary positive/ negative signals and the 

25 derivation of the prior probability in this case is given 
by Hobson & Lasenby (1998) . It is found that the prior 
has the same form as given (3) but the expression for the 
entropy S(s, ro) must be modified. The central idea is to 
express the general (positive/negative) signal vector s 

30 as the difference of two vectors u and v containing only 
strictly non-negative distributions, i.e. 

s=u-v 

By applying continuity constraints on the entropy 
35 functional it is possible that the expression for the 
entropy for the positive/negative signal s is given by 
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S(s, m u , m„) = 2^ - (m u )„ - (m w )„ - s« Jn ^ j j 



(6) 



Where m u and m v are separate models for u and v 
respectively, and where 0n = [s 2 n + 4 (M u ) n (M v )„] 1/2 . We 
cannot hope to replace the models m u and m v by a single 
positive/ negative model m 3 (say) , since such a 
distribution could no longer be considered as an 
integration measure. Nevertheless, we can still consider 
the difference m u - m v as the model for the signal s. We 
note that the form of the positive/negative entropy 
derived by Gull & Skilling (1990) reguires m u = m v . 

Given the entropic prior for general 
positive/negative signals, it is then straightforward to 
define the prior for complex signals simply by applying 
the above analysis to the real and imaginary parts 
separately. In this case the models M u and M v are also 
taken to be complex. The read and imaging parts of m u 
are the models for the positive portions of the real and 
imaginary parts of s respectively, similarly, the real 
and imaginary parts of m v are the models for the negative 
portions of the real and imaginary parts of the image. 
The total entropy of the complex signal is then obtained 
by evaluating the sum (6) using first the real parts and 
then the imaginary parts of s, ra u and m v and adding the 
results. 

The ability to reconstruct positive/negative and 
complex distributions using the MEM approach has profound 
consequences for greatly improving the both the speed and 
accuracy of the MEM technique. These improvements are 
based on the idea of making a change of basis in both the 
signal and data spaces and performing a Bayesian 
reconstruction in this new basis. With an appropriate 
choice of the new bases, it is possible to speed up 
significantly the calculation of the reconstruction, 
which can then easily be rotated back into the original 
basis in signal space to obtain the final reconstruction 
of the signal. For example, the ability to reconstruct 
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complex signals allows us to perform reconstructions in 
the Fourier basis of complex exponentials, which greatly 
reduces the computational complexity of. de-blurring 
images that have been convolved with a spatially varying 
5 point-spread function (see the example below) . 

In order to understand how a general change of basis 
is performed we must first remind ourselves of some basic 
results in linear algebra and vector (see e.g. Riley, 
Hobson & Bence 1997). Suppose there exists a set of 

10 linearly independent vectors e ,n > (n=l,. . . ,N 5 ) , that form 
a complete basis for the N s -dimensional space of the 
signal vector We may then write the signal vector as a 
weighted sum of these vectors. Formally, if we take e (n > 
to be the column vector with unity as the nth element and 

15 zeros elsewhere, then we may write the signal vector as 

N. 



20 



25 



30 



Thus we see that in order to reconstruct the signal 
vector, we are in fact reconstructuring its coefficients 
in this trivial basis. We can, however, equally well 
expand the signal vector in terms of any other linearly 
independent of set of vectors e (nl ' (n=l,...,N s ) such that 



= X>' (n) 



We may perform a similar procedure in the N d - 
dimensional data space, which is in general distinct from 
the signal space. If we consider the trivial basis 
vectors f |n) (n=l,...,N d ) in this space, with unity as the 
nth elements and zeros elsewhere, then the data vector is 
given by Ni 

n = l 

■ 

On performing change of basis in data space to some 
35 other basis f ,n, \ this becomes 

1*4 

n=l 
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Since the noise vector < also belongs to the data 
space a similar change of basis can apply to it, such 
that »„ n< : 



n=l n=l 



It is clear that, even if the element s n of the 
signal vector in the original basis were strictly non- 
negative, the elements s' n in the new basis will in 
general take both positive and negative values. 
Furthermore, in the case where the new basis vectors have 
complex components, the coefficients s' n may themselves 
be complex. Hence it is the extension of the MEM 
technique to the reconstruction of such quantities that 
allows this approach to be taken. 

Once we have performed the changes of basis in the 
data and signal spaces, we denote the vectors with 
components s'„ by s' and we similarly define the vectors 
d« and t' as those containing the elements d' n and < 
respectively. In the signal space we can relate the two 
bases e (nl and e' <n ' (n=l, . . . ,N,)by 

e'00 _ £ Uint V) 

tel 

where U in is the ith component of e' (nl with respect 
to the original (unprimed) basis. The vector s- and s 
are then related by 

s = Us" < 7 ) 
Similar results hold for the bases f (nl and f ' tnl 
(n ■ 1, ,N d ) in data space, such that 

d = Vd' ( 8 ) 
Where the element V in is the with component of f' tn , 
with respect to the unprimed basis. A similar expression 
exists relating the noise vectors c* and <• Substituting 
(7) and (8) and (1), we then obtain 

d' = *'(s') + <'. (9) 
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Where # is a new function relating the signal and 
data vectors in the new basis. 

Clearly, our aim is to choose the new bases in the 
data and signal spaces in order that the relationship (9) 
takes the simplest form. More formally, we wish to 
perform similarity transformations in the data and signal 
spaces that partition each space into numerous (quasi-) 
disjoint subspaces of much lower dimensionality. In such 
bases, we may then calculate the estimate ? of the 
rotated signal vector by performing numerical 
maximisations in each subspace separately. Thus, we 
replace a single maximisation over the N 3 -dimensional 
signal space in the original basis by numerous low- 
dimensionality maximisations over each subspace in the 
new basis. This leads to a considerable increase in the 
speed with which the reconstruction can be performed. 
Then, having calculated the Bayesian reconstruction i' 
in the new basis, the required signal reconstruction s 
can be obtained by rotating back to the original signal- 
space basis, we reiterate that, since the elements S ra in 
the new basis can in general take positive, negative or 
even complex values, it is the extension of the MEM 
technique to the reconstruction of such quantities that 
allows this approach to be taken. 

As an example, let us consider the application of 
the above Bayesian reconstruction technique to the 
deconvolution of images that have been blurred by a 
spatially-invariant point-spread function (PSF) and that 
may also contain some noise contribution. For simplicity, 
we will assume that the de-blurred reconstruction is 
produced on the same pixelisation as the blurred image, 
although this is clearly not required by the technique in 
general. In this simple case, the data and signal spaces 
coincide. 

It is well known that the convolution of an 
underlying image with a spatially invariant PSF is 
equivalent to multiplying together their Fourier 
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transforms and then performing an inverse Fourier 
transform. Therefore, in the Fourier domain, each Fourier 
mode can be considered independently of the others. This 

♦ 

suggests that we should perform the Bayesian 
reconstruction in the Fourier basis, such that 

s' n = J2 exp(-inn{m-\)/N a )s m ^ 

m=l 

With a similar expression relating the components of 
the data vectors d 1 and d and the noise vectors < and t' 
(Since the data and signal spaces coincide) . Thus, in 
this case, the N 3 -dimensional signal (and data) , space has 
been partitioned into N s separate disjoint spaces (i.e. 
one for each Fourier mode) . 

Now for each value of n (or Fourier mode) , we may 
consider the elements d' n , s f n and < independently of 
those for other values of n. This leads to a substantial 
decrease in the CPU time required to de-blur a given 
image. For simplicity, at our chosen Fourier mode we 
denote d' n , s' n and e' n by d', s» and c' respectively. 
The quantity d ! is given simply by the Fourier 
coefficient of the true underlying image, or signal s 1 , 
multiplied by the Fourier coefficient of the PSF, or 
response R. In addition, a noise contribution, c' , in 
the Fourier domain may also be present. If is no 
instrumental errors are expected from a given apparatus, 
it is still possible to introduce "noise" by, for 
example, digitising an image in order to store it on a 
computer. Thus the data value is given by 

d 1 = Rs» + (10) 

Since we are performing the reconstruction in the 
Fourier basis, the noise on each Fourier mode will 
contain contributions from a wide range of scales. 
Therefore, provided the noise on the image is distributed 
in a statstically-homogeneous manner, we would expect 
from the central limit thereon that the noise in the 
Fourier domain is described reasonably well by a Gaussian 
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distribution. Therefore, the likelihood function is given 
by 

Pr(dV) oc exp (-|<f - Rs'\ 2 l° 2 ) . 

where a 7 = (eV) 1 is the variance of the noise 
contribution at the particular Fourier mode under 
consideration. From (6), the entropy S(s» , m) of this 
complex "image" may be shown to be given by (where we 
have set m u = m u = m) 

St.. «) = »(*) - »(«) - ■(') >» [ti^] + * W " 2 ° (m) " 8(0 ln I ' 



where the jr and <p denote the real or imaginary part 
respectively of a complex number; also sr(^) = [S(fi') 2 + ^(m) 2 ! 1 2 

15 and a similar expression exists for O(^). 

Using the above expression for the likelihood and 
prior, and assuming a particular value for the 
regularising parameter d in (5), it is then possible 
numerically to maximise the posterior probability to 

20 obtain the estimator s' of the signal vector at each 

Fourier mode independently. Once these estimators have 
been calculated for all the Fourier modes, we simply 
perform an inverse Fourier transform to recover the 
de-blurred image, i.e. 

25 s n = i- T exp(+iVn(m - lJ/^JC 

The value of a used in the reconstruction algorithm 
may be set in three different ways. Firstly, a may be 
set such that the misfit statistic x 2 > between the 

3 0 observed data and that predicted from the reconstruction, 
is equal to its expectation value, i.e. the number of 
data points to be fitted. This choice is usually 
referred to as "historic" MEM. Alternatively, it is 
possible to determine the appropriate value for o in a 

35 fully Bayesian manner (Skilling 1989) by simply treating 
it as another parameter in our hypothesis space. This is 
the recommended choice and is usually referred to as 
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"classic" MEM- Finally, the simplest option is to fix 

a in advance to some value. This option is unlikely 
to yield optimal results. It can, however, be used to 
obtain a quick solution if the historic or classic value 
5 for * has already been determined for a particular 
problem on a previous occasion. 

The above technique has been applied to several 
different data sets in which an image has been convolved 
with a spatially-invariant point-spread function. For 

10 example, the "classic" de-blurring of a microscope image 
•of a section through a pollen grain. The section has 
dimensions 128 x 128 and is taken from a three 
dimensional "stack" of 20 such images. The original stack 
had been blurred by a spatially invariant three 

15 dimensional PSF and the de-blurred reconstruction of the 
entire stack required approximately 45 seconds on an 
Intel Pentium 200 MH 2 processor. A standard MEM 
algorithm, which does not use the similarity 
transformation technique was also applied to this stack 

20 of images. This produced reconstructions of a similar 
quality to those obtained using the invention, but 
required approximately 50 minutes CPU time on the same 
machine. 

Similar gains in speed can be obtained by expanding 
25 the signal in different bases appropriate to the given 
problem, so long as the correlation matrix of the 
resulting coefficients of the data, signal and noise 
vectors in the new bases is relatively sparse. For 
example, similar results may be obtained by 
30 reconstructing the coefficients in a wavelet expansion 

(Daubechie 1992) of the signal as opposed to the Fourier 
expansion used above. This case has the advantage that 
the coefficients are always real. Furthermore, the 
scaling/ translation properties of the wavelet transform 
35 allow automatic multi-resolution reconstructions of the 
signal vector. 
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Clearly, the general method outlined above can be 
applied to numerous different reconstruction problems of 
arbitrary dimensionality. Examples include the analysis 
of acoustic data, radar, underwater sonar, spectroscopy, 
geophysical data, oil exploration and medical imaging. In 
addition to spatial dimensions, additional dimensions 
such as time, spectral behavior and polarisation are also 
easily accommodated. 
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CLAIMS 

1. A method of reconstructing a signal from a given set 
of data, with a prediction function representing a 

5 predictable effect on the signal, and a noise function 
representing unpredictable noise, the method comprising the 
steps of : 

altering the coordinate basis of the data and signal 
from an original coordinate basis in order to produce a 
10 prediction function having a reduced set of variables; 

performing a Bayesian reconstruction capable of 
operation of positive, negative, and complex signal values 
to produce a reconstruction signal; and 

converting the reconstruction signal back into the 
15 original coordinate basis to generate a signal. 

2. A method according to claim 1, wherein the Bayesian 
reconstruction is performed using a Fourier basis, 

20 3. A method according to claim 1, wherein the Bayesian 
reconstruction is performed using a wavelet basis. 

4. A method according to any preceding claim, wherein the 
Bayesian reconstruction employs the maximum entropy method. 

25 

5. A method according to claim 4, employing an evaluation 
parameter, «, which is determined from a prior 
reconstruction . 

30 6. A method according to claim 4, employing an evaluation 
parameter, oc, which is set at a fixed value. 

7 . A method according to claim 4 , employing an evaluation 
parameter, oc, which is determined during the reconstruction 
3 5 step . 
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8. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is an image signal. 

9. A method according to claim 8, wherein the image 
5 signal is a medical image signal. 

10. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is a radar signal. 

10 11. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is an acoustic data signal. 

12. A method according to claim 11, wherein the acoustic 
data signal is an underwater sonar signal. 

15 

13. A method according to claim 11, wherein the acoustic 
data signal is a geophysical data signal. 

14. A method according to any of claims 1 to 7, in which 
20 the signal to be reconstructed is a signal from 

spectroscopy . 

15. A method according to any one of claims 1 to 7, in 
which the signal is a communication signal, such as a 

25 time-series signal. 
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